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We present a method of filter diagonalization for shell-model calculations. This method is based on the 
Sakurai and Sugiura (SS) method, but extended with help of the shifted complex orthogonal conjugate gradient 
(COCG) method. A salient feature of this method is that it can calculate eigenvalues and eigenstates in a given 
energy interval. We show that this method can be an alternative to the Lanczos method for calculating ground 
and excited states, as well as spectral strength functions. With an application to the M-scheme shell-model 
calculations we demonstrate that several inherent problems in the widely-used Lanczos method can be removed 
or reduced. 



PACS numbers: 21.60.Cs 



I. INTRODUCTION 



To perform numerical investigations of quantum many- 
body systems, many approaches have been proposed, e.g., 
exact diagonalization, the quantum Monte Carlo method, the 
density matrix renormalization group method, and so on. To 
compare with other approaches, the exact diagonalization 
method has a broader range of applications, and can calcu- 
late energies and wave functions without any approximation. 
While a required dimensionality for the Hilbert space is huge, 
the matrix dimension that can be handled in the exact diago- 
nalization approach has recently increased dramatically, ow- 
ing to the development of computers. Hence, the diagonal- 
ization method has become a basic tool in numerical studies, 
and has played an important role in various fields of sciences. 
As for instance, in nuclear structure physics, the exact diag- 
onalization method is of primary importance for shell-model 
calculations. 

For an exact diagonalization in large-scale shell-model cal- 
culations, the Lanczos method [1] has so far been the only 
feasible method for practical use. This method has been 
widely employed to obtain not only ground states but also 
low-lying excited states. Nevertheless, there still exist three 
long-standing problems: (1) In calculating highly excited 
states, convergence is much slower than that for the ground 
and low-lying states. The number of the Lanczos iteration 
process tends to grow rapidly as the energy goes higher. (2) 
The Lanczos method needs to do reorthogonalization of all 
obtained Lanczos vectors, which demands substantial numer- 
ical effort. This problem is rather technical but crucial in prac- 
tice because the reorthogonalization procedure sets a practical 
limitation in solving highly excited states. (3) In large-scale 
shell-model calculations with the M-scheme, the total angular 
momentum J and the total isospin T are not necessarily con- 
served for each basis, although the total magnetic quantum 
number J z = M is conserved by definition. Then conservation 



of angular momentum and isospin may be violated in some 
cases. In the Lanczos method, conservations of J and T can 
be realized by choosing an initial wave function with good 
quantum numbers J and T, However, this procedure is not so 
stable against round-off errors. Therefore, the conservation of 
these quantum numbers is an important issue particularly in 
the M-scheme shell-model calculations. 



Up to now, several shell-model codes [2-4] have been de- 
veloped for state-of-the-art large-scale calculations. However, 
there has been no attempt to solve the long-standing and basic 
problems in the Lanczos method mentioned above. 

Recently Sakurai and Sugiura (SS) JH |(| have proposed 
a new diagonalization method for a generalized eigenvalue 
problem: Ax = XBx, where A and B are arbitrary matrices 
(i.e., not necessarily symmetric matrices). Their method is ap- 
plicable even to complex matrices. In this method, Cauchy's 
integral formula is used in order to obtain eigenvalues (and as- 
sociated eigenvectors) inside of the region enclosed by a given 
integration contour, which can be considered to be a kind of 
a filter. Therefore we call this new method A_gfilter diagonal- 
izationA_h hereafter. 

In the SS method, a diagonalization problem turns into a 
problem of solving a large number of linear equations, which 
also demands a heavy computation for large-scale shell-model 
calculations. To overcome this difficulty, we use the shifted 
complex orthogonal conjugate gradient (COCG) method [7]. 
The shifted COCG method corresponds to a combination of 
"shift" algorithms H] and the COCG method @], which is 
designed to solve a particular family of linear equations. An 
advantage of the shifted COCG method is that a problem of 
diagonalization can be reduced to just one linear equations. 
With the help of the shifted COCG method, the SS method is 
greatly reinforced and becomes more feasible. The first study 
on the SS method with the shift algorithms was presented in 
Ref. fioll . Very recently, an application and an extension of the 
SS method with the shift algorithms have been reported for 
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all-to-all propagators in the lattice quantum chromodynamics 

(Qcd)Gj. 

In this paper, we apply the filter diagonalization based on 
the SS method combined with the shifted COCG to quantum 
many-body systems, and demonstrate that the filter diagonal- 
ization is indeed an alternative to the Lanczos method in eval- 
uating energy eigenvalues, eigenstates and spectral strength 
functions. Moreover, the aforementioned problems of the 
Lanczos method in the M-scheme shell-model calculations are 
shown to be removed or reduced. 

This paper is organized as follows: In Sec. II, we show 
the filter diagonalization based on the SS method and the 
shifted COCG method, and present how to evaluate the spec- 
tral strength function. In Sec. Ill, we present several examples 
of numerical calculations and discuss characteristic properties 
of the method. In Sec. IV, we give a conclusion. In Appen- 
dices, we summarize useful relations concerning the Hankel 
matrix and an algorithm of the shifted COCG method. For 
readers who have interest in this diagonalization, this paper is 
written in a self-contained manner. 



II. FILTER DIAGONALIZATION OF SHELL-MODEL 
CALCULATIONS 

A. SS method 

In this section, we summarize the SS method in the shell- 
model calculations. In order to reduce a large-scale eigen- 
value problem to a small scale one, we first consider moments 
fl p {p = 0, 1,2, • • • ) defined by Cauchy's integral as, 



1 r (z-e) p 



(II 



where \\jr) and |0) are arbitrary wave functions, and H is a 
shell-model Hamiltonian, satisfying the eigenvalue equation 
H\(pi) = ei\<Pi). £ denotes the energy in the vicinity of an en- 
ergy region of interest (target region). F means an integration 
contour to enclose energy eigenvalues in the target region, as 
depicted in Fig.[T] The integration is carried out on the com- 
plex z plane, so that energy eigenvalues on the real axis are 
energy poles if they are inside the integration contour F. As 
a result, these eigenvalues contribute to the integral, and they 
are central quantities in the SS method |5|]. 
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FIG. 1: (Color online) An illustration of integration contour F (an 
open circle) and energy poles (filled circles) on the complex z-plane. 
In this illustration, F encloses one of the energy poles on the real-z 
axis. 



To clarify the physical meaning of these moments, we ex- 
pand \y) and \(j>) in terms of the ortho-normalized energy 
eigen functions \(p)'s of the Hamiltonian H, that is, = 
^Cj\(pi) and |0) = ^di\(pi), where c's and a"s are coefficients 

with £| Ci | 2 = land £|<4| 2 = 1. 

Due to the theorem of residue, Cauchy's integral is formally 
carried out and the moments are rewritten as 



Hp 



e) p c k d k . 



(2) 



The summation over k is taken if energy eigenvalues are inside 
the r. The moment jl p vanishes when none of the energy poles 
is enclosed by F, or when amplitude is zero for the eigenstates 
corresponding to the poles (i.e., c^dk — 0). 

To extract the energy eigenvalues e k (k e F) from these mo- 
ments, we follow the SS method 

H. Namely, we solve the 
generalized eigenvalue problem formulated as 



Mx = XNx, 

where M and N are the n x n Hankel matrices defined by 



(3) 



M = 



/Mi. M2, 
' Hz, M3, 
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N = 



( Ikh Mb 
ill, flz, 



V Un- 



it V-n, 



M2n-1 / 



Mn-1 \ 



M2n-2 / 



(4) 



(5) 



It is then possible to demonstrate that the eigenvalues X k in 
the generalized eigenvalue equation correspond to e k — e. Its 
proof needs a property of the Hankel matrices that they can be 
always factorized with the Vandermonde matrix as shown 
in Appendix A. Note that this method to extract eigenvalues 
from moments was used in Ref. II 1311 . 

The dimension n introduced in the generalized eigenvalue 
equation corresponds to the number of eigenvalues inside the 
integration contour, but it is not known a priori. The optimum 
n can be obtained by monitoring a convergence pattern of the 
energy eigenvalues as a function of n. This is because the 
energy eigenvalues should be unchanged when the n exceeds 
the number of eigenvalues inside the integration contour. 

The amplitude c k d k of (e k — e) p in Eq. (f2]i can be obtained 
by the diagonal matrix given as 

D = V- 1 N(V T )-\ (6) 

where V is a Vandermonde matrix defined by Vij = (ej — 
e)'- 1 , i.e., 

\n-l 



( 1, ei-e, ■•• (ei-e)"- 1 \ 
1, e 2 -£, ■•• (e 2 -£)"- 1 



V 1, e n -e, ■■■ (e n -e)"- 1 J 



(7) 
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because of Eq. (IA6l l in Appendix A. It should be noted here 
that inverse operations of the Vandermonde matrix are not nu- 
merically stable. It is thus better to use the eigenvectors in 
the generalized eigenvalue equation in practical calculations, 
because (V T )~ l is equivalent to the eigenvectors of Eq. (0. 

To study electromagnetic transition properties, wave func- 
tions should be described in the framework of the SS method. 
For this purpose, we define vectors \s p ) as 



2niJ T z—H IT, "~' (8) 
In the same way as Eq. (O, this can be also formally rewritten 



\s p )=Y,M<h){V T )k P - (9) 
keT 

Therefore, wave functions |^) are explicitly obtained as, 



l^)-LW(v r ) 



T\-\ 

pk ■ 



(10) 



Its general proof is shown in Ref. 

An error analysis was presented for the Hankel and Van- 
dermonde matrices in the context of the SS method, in Refs. 
C3. 



B. Numerical integration, scaling, and shifted COCG method 

Next, we explain how to integrate the moments jj. p . An 
integration contour F is chosen to be a circle given as, 

z = e + re ie (e,r:real, 6 = [Q,2n]). (11) 

The target eigen energies are then located between e — r and 
e + r. Cauchy's integral is now evaluated numerically by the 
trapezoidal rule with respect to angle 9 as 



Aip^ir L v — — b— *>> 

No k =o z k -H 



(12) 



where zt = £ + re 7 ^ + - . Here, we take integral points in 
a symmetric manner about the real axis because we take an 
advantage of the property f(z) = f(z) for a complex number 
Z. 

The integration contour with a larger r can include more 
energy poles. However, as the ji p has the r p dependence, the 
moments become larger as a function of p, which causes a 
numerical instability in Eq. (01. To remove it, we scale Eq. (Q3 
by mapping the circle with radius r into a unit circle Jfl as 



z' = z/r = e/r + e' 6 . 
Under this mapping, the moments \l' become 

ker r 



(13) 



(14) 



where jl' p = \l p jr l '■ Then, the r p dependence is removed in 



1 



N -l 



iV ° k=0 



(4-0 



(15) 



where z' k — Zk/r, H' — H/r and e' = e/r. 

For each angle 0, we need to evaluate a matrix element 

(v| — which involves an inverse operator. To avoid 

handling inverse operators, we define \%) as 



\<j>) = (z-H)\ X ), 
and calculate \%) first, then obtain (y/\ — 



(16) 



H 



10} = (viz)- 



To obtain \%), we solve linear equations; Ax — b, where 
Kn.n = (m\z—H\n), b m = (m\<j>) and x m = {m\%). A vector 
to) means an M-scheme basis. This equation is solved by the 
COCG method |0] for complex, symmetric, but non-hermitian 
matrices, because complex number z appears in the diagonal 
matrix elements. As \%) depends on z, the above linear equa- 
tions should be solved for each z- As the number of integral 
points A^o increases, this numerical calculation becomes more 
time-consuming. However, by using an invariance property 
of the Rrylov subspace, we can drastically reduce the amount 
of computation. Once we can solve \<j>) = (zq —H)\xo) at a 
certain zo by the COCG method and store residual vectors, we 
can compute 1 } = (z — H) \ % ) for z ~ zo from the stored resid- 
ual vectors. This method is called the shifted COCG method 
10, Q3]. Details are shown in Appendix B. We will present 
how to reduce computation by this method in Sec. III. B. 



C. Spectral strength function 

To investigate a dynamic property of a system concerning 
an operator O, it is useful to evaluate a spectral strength func- 
tion I(co) defined as, 

/M = Ll(v/i fi Vl^ ) }| 2 5(c-(£i fi) -4 A) )), d7) 

n 

where and are energies of the n-th state and the 0- 

th state, respectively, and ly^ 5 ') and ly^') are the associated 
eigenstates. If the operator O violates the conservation of cer- 
tain quantum numbers, e.g., angular momentum, isospin, and 
numbers of proton and neutron, the initial and the final states 
can belong to different Hilbert spaces indicated with labels 
A and B. By a relation 1/ (x+ ii]) = P[l/x] - in8(x), the 
strength function can be rewritten as, 



, , 1 

I (CO) = Im 

% 



1 



CO 



,(A) 



H + irj 



■o\% 



(Ah 



(18) 

where rj means a half width. Here we define a complex num- 

(A) 

ber z as z — CO + E^ +irj and a new normalized wave function 
belonging to the B space as, 



\^) = o\tf ) )/y/(tf ) \oio\vi? ) ). 



(19) 
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Then, evaluation of the strength function can be reduced to 

(B) 1 (B) 

the calculation of the matrix element (<Pq | Jj^o )• By 

the Lanczos method, the Hamiltonian matrix is transformed 
into a tridiagonal form with matrix elements which are usually 

(B) 1 (B) 

denoted as a,- and /3,. The matrix element (<Pq | \(p^ ) 

can be expanded in the form of continued fraction ifisjf as 



In practical applications, as various properties of wave func- 
tions are also important, we often calculate the eigenstates in 
addition to the eigen energies. In such cases we can directly 
evaluate the strengths by using Eq. d!71 l. which is equivalent to 
Eq. ( f20b - The half width is also introduced by the Lorentzian 
curve. 

{ B\ 

By the Lanczos method starting from \ (p^ ), strength func- 
tions converge faster as z becomes smaller. To obtain the 
strength function of higher excitation energy, the number of 
the Lanczos iteration is increased inevitably, which results in 
a serious "inflation" of computation time for matrix elements 
calculations and the I/O access time to storage devices due 
to the reorthogonalization among the Lanczos vectors. More- 
over, in the M-scheme calculations for large-scale shell model, 
the Lanczos method often fails to conserve angular momen- 
tum through numerical errors, so that a delicate treatment is 
necessary for their conservation as will be discussed later. In 
general, such calculations are quite difficult. 

Next we consider the filter diagonalization for the spectral 

strength function. To obtain excitation energies E n B ^ — 

and matrix elements {y^\0\ y^}, two states \ \\r) and \<j>) in 

Eq. ([T]) are set to be <9|y/^ A '). By expanding 0| V^ A ' ) with the 

complete set m the B space as 0|v^ A ') = 52 N V 7 ,^)* 

the moments in Eq. (O are rewritten as 

H p = Y,{E { n ] -e) p b 2 n (21) 

where b 2 = \{Wn B \o\^f^)\ 2 . By the filter diagonalization, 

we can obtain E n B ^ and b 2 due to Eq. ©, and therefore we 

can plot b 2 as a function of excitation energies E„ B ^ — E^\ 
Compared to the Lanczos method, it is advantageous that we 
can directly evaluate the strength function in a given excita- 
tion energy region. Moreover, aforementioned problems in 
the Lanczos method are removed or reduced, which is demon- 
strated in Sec. III. E. 

III. NUMERICAL TESTS 

A. Lanczos method and conservation of quantum numbers 

To test the filter diagonalization in the shell-model calcula- 
tion, we consider 48 Cr in the model space consisting of single- 



particle orbits /7/2)/?3/2 >/5/2 an d Pi/2- Its M-scheme dimen- 
sion for M = is about 2 x 10 6 . This calculation used to be 
a state-of-the-art large-scale shell-model calculation in 1994 
ifjfSll . so that it has been often used as a benchmark test for 
new shell-model methods 1 17-19]. Moreover, due to N = Z, 
the M = space contains all states with angular momentum 
0,1,2, - • • and isospin 0,1,2, - • •. It is a touchstone whether the 
filter diagonalization can handle such quantum numbers cor- 
rectly. In this work, we use the KB3 interaction rf20Tl as a 
residual interaction. 

In the large-scale shell-model calculations, the M-scheme 
is often used but it has a problem in conservation of angular 
momentum and isospin. In principle, conservations of / and T 
should be maintained if we take an initial state with good J and 
T, but it works well only for simple cases. For instance, let 
us suppose the Lanczos iteration, starting from an initial state 
with / = 0. It is easy to obtain a ground-state wave function 
having 7 = 0, but it is not so for excited states. This is because 
numerical round-off errors can give rise to eigenstates with 
different angular momentum. 

For such a case, we can manage to deal with this problem by 
introducing a modified Hamiltonian H' = H + aJ ■ J + fiT ■ T 
with positive a and j3, which push up undesired components 
into higher energy region. Although this technique is widely 
used and works well, it is applicable only to ground and low- 
lying states. 

Higher excited states with J — are quite difficult to obtain 
by the above approach, because the M — space also contains 
states with non-zero angular momentum 7^0. Small numer- 
ical round-off errors can easily contaminate the 7 = wave 
function with wrong components (7 ^ 0). In such a case, the 
double Lanczos method 112 ill has been proposed. That is, in 
addition to the usual Lanczos iterations for each Lanczos vec- 
tor, we apply the Lanczos diagonalization concerning the 7 • 7 
term (and T ■ T). This additional Lanczos process can remove 
the unnecessary components of non-zero angular momentum 
(and isospin) caused by the round-off errors. 

In Figs. 3 and 4, we show the lowest 12 energies of 7 = 
and T = states calculated by the double Lanczos method. 
Table I is a list of the numbers obtained by the two kinds of 
iterations. The number of the main Lanczos iterations and 
the total number of additional Lanczos iterations for J ■ J are 
denoted as Nl(H) and Nl{J 2 ), respectively. For the excited 
states with 7 = 0, the double Lanczos method starts from the 
lowest 7 = state in the {/t/2) & configuration space. For the 
ground state, Nl{J 2 ) is zero as expected, while Nl(J 2 ) rapidly 
increases for higher excited states. In this way, it was demon- 
strated here that the double Lanczos calculation needs addi- 
tional (and heavy) computational efforts. Nevertheless, there 
had not been a better way than the double Lanczos method, 
so that it was inevitably an indispensable approach in obtain- 
ing excited states with good J in the M-scheme shell-model 
calculations. 
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N L (H) 


17 


30 


38 


47 


163 


Nl(J 2 ) 





17 


50 


88 


668 



TABLE I: The number of the main Lanczos iterations and the to- 
tal number of additional Lanczos iterations for J • J are denoted as 
Ni(H) and Ni(J 2 ), respectively. They are calculated for several low- 
est states with J = and T = of 48 Cr. 



B. Test of ground and low-lying states by the filter 
diagonalization 

Next we consider the filter diagonalization in the shell- 
model calculations. 

First of all, we calculate the yrast states of 48 Cr at J = 0,2,4 
and 6 as an example, with an aim to demonstrate how the filter 
diagonalization is proceeded numerically. To evaluate the mo- 
ments defined by Eq. (Q]), arbitrary states \<p) and | y/) need to 
be prepared. In the original SS method, they were chosen to be 
vectors consisting of random numbers. Instead, here we em- 
ploy lowest energy wave functions obtained through a diago- 
nalization of the Hamiltonian matrix in the two-particle two- 
hole (2p2h) space, i.e., {fyif^iPiliJsiiiPi/iY < 2). 
These wave functions are approximated states with good an- 
gular momentum (J = 0,2,4,6) andisospin (T = 0). Hereafter 
we call these states for \<p) and | y/) "initial states" in the con- 
text of the filter diagonalization. The dimension of the M = 
(2p2h) space is 62220, while the dimensions of the M + 1 
spaces are smaller. These cases can be easily solved by means 
of the standard diagonalization techniques. The energy of the 
lowest state with J = is -31.1 MeV. 

As for an integration contour F, we take a circle with radius 
r, which covers an energy interval [e — r, £ + r] . In Fig. [2] we 
choose a different circular integration contour for each J, of 
which center is at z = -33.0 MeV for 7 = 0, -32.0 MeV for 
J = 2, -3 1 .0 MeV for J = 4 and -29.5 MeV for J = 6. The 
radius r is 1.0 MeV. These integration contours cover energy 
intervals [-34, -32], [-33, -31], [-32, -30] and [-30.5, 
—28.5] (in MeV) for/ = 0,2,4 and 6, respectively. Numerical 
integrations are carried out by means of the trapezoidal rule. 
As shown in Fig. [2] ten points along the contour are used for 
the numerical integration. ( Note that in practice it is sufficient 
to calculate only at five points located in Imag(z) > 0, due to 
the property /(z) =f(z).) 

For numerical evaluations of the moments, at each point 
on the integral contours, it is possible to solve a set of linear 
equations, Eq. ( fTBT ) by means of the COCG method. This cal- 
culation, however, tends to be quite time-consuming, as the 
number of integral points increases. To reduce the amount 
of computation, we use the shifted COCG method. With the 
shifted COCG method, once we solve \(j>) = (zq —H)\xo) for 
a particular zq, solutions at the other neighboring points z ~ zo 
can be obtained with a small computational cost, if the itera- 
tion number needed for the convergence at z is less than that 
at zo. This condition will be discussed later. First, Eq. (fT6] > is 
solved atzo = -32.0 + 0.1/ MeV for the J = state. The solu- 
tion was obtained by 19 iterations under the convergence crite- 
rion that the norm of the residual vector is less than 10~ 5 . The 



values of the integration at other integral points for the / = 
state are obtained by the shifted COCG method. Therefore 
the computational cost does not nearly depend on the number 
of integral points. It mainly depends on the iteration numbers 
of the COCG method at zo- Thus exact ground state energy 
is obtained by this filter diagonalization with almost the same 
computational cost as that of the Lanczos method (see Table 
I). 

In Fig. 12 the integration contour To encloses two energy 
poles for 0i and 2j states because of the M = space. How- 
ever, as we always use an initial state with good J, eigen states 
with different J can be filtered out and such states never ap- 
pear in the solutions of Eq. (f3j). 
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FIG. 2: (Color online) Demonstration of the filter diagonalization 
for the yrast states of 48 Cr on the complex z-plane. The yrast- 
state energies obtained by the filter diagonalization and the Lanczos 
method are shown by crosses and small circles, respectively. For 
J = 0,2,4 and 6, the COCG method is applied at z = -32.0 + 0.1/, 
-31.0 + 0.1/, -30.0 + 0.1/ and -28.5 + 0.1/ (in MeV), respectively 
(diamonds). The numerical integrations are carried out separately for 
each angular momentum using 10 points along the contour, which are 
shown by squares. Horizontal and vertical axes are real and imagi- 
nary parts of z, respectively. 

Next we consider the low-lying excited states with 7 = 
and T = quantum numbers. In Fig.[3](a) and (b), circles with 
r = 1 and r = 2 are shown, respectively, which cover the same 
energy interval [—33.5, —24.0] in MeV. In these calculations, 
we take the lowest state in 2p2h space as an initial state in 
Eq. ©. 

Figure[3](a) is an extension of Fig. [2] for the 7 = state in 
wider energy regions. Numerical integration is carried out by 
20 points for each circle, and we carry out the COCG calcu- 
lation only at z = 24.5 + 0.1/. For the other integral points, 
the values of the integrand are obtained by the shifted COCG 
method. For the circle with a center at z = —30.7 MeV, the 
moments vanish. It means no eigenvalue in this energy inter- 
val [—31.7, —29.7] in MeV. In the following circles, we can 
confirm the energies for 0i, 02,03, and O5 states. Because the 
initial state is / = and T = and matrix-vector multiplica- 
tions in the COCG method conserve the quantum numbers, 
no state with different quantum numbers appears. Compared 
to the Lanczos method, the filter diagonalization is found to 
be advantageous with respect to the conservation of quantum 
numbers in numerical calculations. 
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FIG. 3: (Color online) Demonstration of the filter diagonalization 
for excited states of 48 Cr on the complex z-plane. The convention of 
symbols (crosses, circles, diamond and squares) is the same as that of 
Fig.0 The energies located in [—33.5, —24.0] MeV are calculated 
using two integration contours with different radii, (a) r = 1 MeV and 
(b) r = 2 MeV. Horizontal and vertical axes are real and imaginary 
parts of z, respectively. 



In Fig. [3(b), we use circles with radius r = 2 MeV, which 
give us the same results. In this calculation, we use Eq. (15[ 
for scaling. For both calculations, O4 state is not reproduced 
because the initial state has very small components of the O4 
state (0.03%). 

In Fig. [4] energy interval [—27.5, —22.5] in MeV is shown. 
Here we use smaller circles with r = 0.5 MeV. Because a 
smaller circle includes fewer eigenvalues, it is easier to solve 
the equation. Smaller circles are expected to be useful when 
the level density is large. However, as shown in the next sub- 
section, convergence of the COCG method unfortunately be- 
comes slower. 

In this calculation, as an initial state, we use the sum of 
the lowest five wave functions with / = and T = in the 
2p2h space and can reproduce 03^9 states, including the O4 
state. As shown in Eq. (O, since Cauchy's integral makes use 
of an initial state to extract eigenstates within the integration 
contour, the choice of the initial state is important. 



C. Convergence of the COCG method 

The computational cost of the present method mainly de- 
pends on the convergence property of the COCG method. Due 
to the shifted COCG method, dependency on the number of 
integral points or the size of the integration contour is very 
weak. In Fig. [5] we show several convergence patterns of the 
COCG method at z= -31+0.1/, -27 + 0.1/ and -23 + 0.1/ 
(in MeV). Here the norm \r\ of the residual vector defined 
in Eq. ( IB4l i is plotted as a function of the number of itera- 



3 



□ 


□ □ □ □ □ 


° □ n D 


1 

° □ □ □ 


° □ n D 


° □ □ n 


d □ 


□ 


□□ 


□□ 


□□ 


□□ 


□□ 


nd 


□ 


□□ 


□□ 


□□ 


□□ 


□□ 


□□ 




□ □ 


□ □ 




® J 9» 








» □ 3» 


□ □ 


□ □ 




□ □ 




□ 


□ □ 


□□ 




□□ 


□□ 


□□ 


□ 


□ □ 


□□ 


□ □ 


□□ 


□□ 


r. na 


□ 


□ □ □ □ □ 


□ □ □ □ 


□ □ □ □ 

1 


□ □ □ □ 


□ □ □ □ 


□ 



-27.5 



-25 
Realz (MeV) 



-22.5 



FIG. 4: (Color online) Demonstration of the filter diagonalization 
for excited states of 48 Cr on the complex z-plane. The convention of 
symbols (crosses, circles, diamond and squares) is the same as that 
of Fig. [2] The energies located in [—27.5, —22.5] MeV are solved 
using integration contours with r = 0.5 MeV. Horizontal and vertical 
axes are real and imaginary parts of z, respectively. 



tion of the COCG method. We take |r|/|fe| < 10~ 5 as a crite- 
rion of convergence, and as an initial state, we take the sum 
of the lowest five wave functions with J = and T = in 
the 2p2h space. In general, the convergence pattern of the 
COCG method is not monotonic, but on average, the norm of 
the residual vector decreases. As real part of z increases, the 
number of iteration for convergence increases. 

To investigate the z dependence of the number of iteration, 
in Fig. [6j its contour plot on the complex z-plane is shown. 
The energy eigenvalues are also shown on the real axis by 
open circles. In general, as imaginary part of z increases, the 
number of iteration decreases. As real part of z increases, 
the number of iteration also increases. Along a given inte- 
gration contour, the number of iteration of the COCG method 
becomes largest at the point z whose real part is largest and 
imaginary part is smallest. Therefore, in Figs. [2HH such a 
point is chosen as the zo of the COCG method, and the values 
at the other integral points are obtained by the shifted COCG. 

Globally the COCG method converges fast for the ground 
and several low-lying states, while its convergence becomes 
worse for highly excited states. For such energy eigenvalues 
further theoretical development is necessary. 



D. Numerical accuracy 

In the filter diagonalization, we use numerical integration 
to evaluate the energies and wave functions. Here we discuss 
their numerical accuracy. For example, we again consider 
the calculation of the ground state, taking the lowest state in 
the 2p2h space as an initial state. The ground-state energy is 
—32.954 MeV. Like Fig. [2j we enclose this energy pole by 
one circle with radius r = 1 .0 MeV. By moving its center po- 
sition e, the ground-state energy pole is located at center or 
peripheral of the circle. 

In Fig.|7](a) and (b), we plot the moments /Xo, Hi and energy 
as a function of the center position e for two cases of 10 and 30 
integral points. Here flo is a square of an overlap between the 
initial 2p2h wave function and the ground state, and energy is 
given by the ratio of these two moments, /jiio + £ because 
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FIG. 5: (Color online) Convergence patterns of the COCG method. 
The norm of residual vector is shown as a function of iteration num- 
ber for z = -31 +0.1/, -27 + 0.1 i and -23+O.li (in MeV). 
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FIG. 6: (Color online) Contour plot of the iteration number of the 
COCG method on the complex z-plane. Horizontal and vertical axes 
correspond to real and imaginary part of z, respectively. The energy 
eigenvalues are shown by open circles on the real axis. 



the integration contour encloses one energy pole. In Fig.[7](b), 
the energy is quite constant as a function of the center position, 
although at e = -32.954 ± 1 .0 MeV, the moments should be 
divergent, and for e < -33.954 MeV or e > -31.954 MeV, 
both Hq and ji\ should be zero. 

On the other hand, in Fig. 0(a), we can see that each mo- 
ment ill-behaves at such critical values. From Eq. (0, jUn is 
constant and pi\ = (—32.954 — e)jJ^. When the energy pole 
comes to the peripheral of the circle, jio deviates from a con- 
stant value and /Xi does not follow the linear behavior. By 
increasing the number of integral points, we can see that the 
numerical accuracy is improved. However, when the obtained 
energy is close to e ± r, the energy itself may be still valid but 
the absolute values of the moments lose their reliability. 

Next we consider the reliability of the calculation of wave 
functions. By using Eqs. (O and ( [Tol l, we can explicitly obtain 
wave functions. In Fig.[7](c), we plot the overlap between the 
ground state wave functions obtained by the Lanczos method 
and by the filter diagonalization as a function of the center po- 
sition. The overlap is also quite constant like energy. The ill- 



behavior comes from the denominator which can change the 
norm of wave functions. By renormalizing the wave function, 
this ill-behavior can be weakened. In Fig.[7](c), the overlap 
between the initial state and the ground state obtained by the 
filter diagonalization is also quite constant. As this quantity 
is the same as ;iio, the /J.q obtained from the wave function is 
more reliable. Thus in the filter diagonalization, the accuracy 
of energy and wave function is better than that of the absolute 
values of the moments. 



(22) 



Note that, by the energy variance 0~ [ 1 9] defined as 



(H 2 )-(H) 2 



(H) 2 



we can evaluate the quality of the calculations without any ref- 
erences. In this case, this a is perfectly zero, which means that 
the obtained energy and wave function are exact. The compu- 
tational cost of a is the same as that of the energy expectation 
value and this a can be easily numerically evaluated. 
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FIG. 7: (Color online) The moments (a), energies (b) and overlaps (c) 
are plotted as a function of the center energy of integration contour, 
which is a circle with radius r = 1 .0 MeV. Two results for 10 and 30 
integral points are shown by dotted lines with open circles and solid 
line with filled circles, respectively. In (b) and (c), two results are 
almost the same. In (c), upper line with marks (sky blue) shows over- 
laps between the ground states obtained by the Lanczos method and 
by the filter diagonalization. Lower line with marks (green) shows 
the same quantity as the but it is evaluated by the obtained wave 
functions. 
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E. Test of Ml strength function 

As an accuracy test for the spectral strength functions ob- 
tained with the filter diagonalization, we consider an Ml 
strength function of 48 Cr. 

The ground state \\j/o) with 7 = and T = is obtained 
by the Lanczos method or the filter diagonalization. The Ml 
operator O with the free g-factors is given as 

O = g?L* + gjL v + g*S* + g v s S v (23) 

where L K and L v are the proton and neutron orbital angu- 
lar momentum operators and S K and S v are the proton and 
neutron spin operators, respectively. The free g-factors are 
gf = 1, gj = 0, gf = 5.586 and g* = -3.826. We consider 
the |<po) = 0\Yo), of which angular momentum is 1, while 
the Ml operator O mixes isospin. Then we classify the Ml 
operators as, 

O = O T=0 + O T=1 (24) 

and 

t=o = gj Igj (L"+I?) + 8s 2 gs (S K + S V ) (25) 
Q T=i = Si ~8i ( L x_ L Vj + 8s-8s ( S n_ s vy (26) 

As an initial wave function of the filter diagonalization, we 
prepare | <pt>) = T=0 \yo) and <p ) = <9 r=1 |v/ ), of which 
angular momentum are 1 and isospin are and 1, respec- 
tively. By this technique, the filter diagonalization is carried 
out within the specified space. 

In Figs.|8](a)-(c), we present several strength functions ob- 
tained by the double Lanczos method with different numbers 
of Lanczos iterations. Lower energy part of the strength func- 
tion converges fast as a function of the number of Lanczos iter- 
ations, while convergence of higher energy part of the strength 
function is slow. In Fig.[8](d), we present the results of the 
filter diagonalization. We can see that the present filter diag- 
onalization can correctly reproduce the Ml strength function, 
compared to Fig.|8](c). 



IV. CONCLUSION 

In this paper, based on the SS + shifted COCG method, we 
have shown an alternative diagonalization method for shell- 
model calculations. This method is called the filter diagonal- 
ization. It has a salient feature that eigenvalues and eigenstates 
can be searched for within a given energy interval. The filter 
diagonalization works equally well or is superior to the Lanc- 
zos method. Since both methods are based on the property of 
the Krylov space defined by eq. dBl U . their basic frameworks 
are similar. However, the following differences can distin- 
guish one from the other. 

In state-of-the-art large-scale shell-model calculations, the 
M-scheme is very useful but it needs a delicate treatment for 
angular momentum and isospin. In the numerical calculations, 




Excitation energy (MeV) 



FIG. 8: (Color online) The Ml strengths divided by the total strength 
as a function of excitation energy. The results are obtained by the 
double Lanczos method with (a) 50, (b) 100, (c) 500 iterations, while 
(d) by the filter diagonalization. The curves show the results of fits 
by a Lorentzian with a half width of 200 keV. 



the robustness of conservation of these quantum numbers is 
different between the two methods. In the Lanczos method, 
small round-off errors easily break down such conservation, 
so that the double Lanczos method [21] was developed. On 
the other hand, in the filter diagonalization, conservation of 
the quantum numbers is found to be quite robust, which is 
a superior property. One of the problems in the Lanczos 
method, when applied to large-scale calculations, is reorthog- 
onalization of the Lanczos vectors, which demands a heavy 
I/O access to storage devices. In the filter diagonalization, we 
use residual vectors, which are similar to the Lanczos vec- 
tors, but reorthogonalization is not necessary. This is another 
superior property. Because of the two merits, the filter diago- 
nalization is superior to the Lanczos method especially for the 
calculations of excited states and spectral strength functions. 

To examine such properties of the filter diagonalization, we 
have investigated its feasibility by taking 48 Cr as an example 
with the configuration space consisting of / 7 pi, p-} pi, fs pi, and 
Pin orbits. This calculation is often considered as a touch- 
stone of a new method aiming at large-scale shell-model cal- 
culations. We have demonstrated that while keeping good an- 
gular momentum and isospin, the filter diagonalization can 
obtain the yrast states and off-yrast states efficiently, and that 
it can also be useful for spectral strength functions. As for 
larger-scale calculations, we have tested the filter diagonal- 
ization for the case of 56 Ni with GXPF1A interaction M22I1 . 
The 8p8h space [23] has approximately 2.5 x 10 8 dimension. 
We can correctly obtain the ground state, oblate and prolate 
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deformed states by the filter diagonalization. 

Finally, we point out two open problems. One is the con- 
vergence of the COCG method, which depends on the position 
of complex energy z. For highly excited states, convergence 
becomes slow. The other is how to choose the integral con- 
tour and integral points for more efficient or unskilled com- 
putation. The present integral contour is circle but this is not 
unique [11]. Other integral contours may be more convenient 
and may solve the convergence problem. For these problems, 
further theoretical developments are strongly needed. 
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where 



A = 



M = VDAV 1 , 



/a h 0, ••• \ 
0, a 2 , ••• 

V 0, 0, ••• a n I 



By these factorizations, we can prove JH] 



(A7) 



(A8) 
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M - XN = VD(A - Xl)V 7 



(A9) 



Therefore, eigenvalues of generalized eigenvalue equation, 
Mx = XNx, are A = a k (k = 1,2,3, ••• ). 



Appendix A: Factorization of Hankel matrix 

Here we summarize a factorization of the Hankel matrix. 
The moments are defined as 



(Al) 



where a k and b k are, in general, complex numbers. The n xn 
Hankel matrix is defined as 



N = 



/ jU , Ml, ••• Mn-l \ 
Mb M2) ••• Mn 



V Mn-b M") ' ' ' M2n-2 / 

I Lb k , lakh, ■■■ ,Z<4r 1 h \ 
\b k , ■■■ ,lLa n k 



JM k b k , T.a 2 k b k , ■■■ Xa" k b k 



(A2) 



(A3) 



\Za n k - l b k: La" k b k , •■■ M k b k J 

The n x n Vandermonde matrix V and diagonal matrix D are 
defined as 



V T = 



and 



D = 



I 1, a u 

1) «2, 



\ 1, a„, 



fb u 0, 
' 0, b 2 , 



» 



(A4) 



(A5) 



V 0, 0, ••• bj 
Therefore, the following factorization holds as, 
N = VDV T . 



(A6) 



Appendix B: Shifted COCG method 

The conjugate gradient (CG) method is an algorithm to nu- 
merically solve linear system as 



Ax - 



(Bl) 



where A is a matrix and x and b are vectors. We consider the 
following quadratic function f(x) defined as 



f(x) = ^x T Ax — x T b. 



(B2) 



At the stationary point x m , where f'(x m ) — 0, the equation 
Ax m = b is satisfied. Therefore, we iteratively minimize f(x) 
by changing x along negative gradient direction, starting from 
Xq. A merit of the CG method is that we can handle only 
multiplication of matrix A to vector x. During iteration pro- 
cess, matrix A is unchanged and sparseness of matrix A al- 
ways holds. In the application of quantum systems, it is very 
useful for conservation of quantum numbers. 

The complex orthogonal conjugate gradient (COCG) 
method [9] is a generalization of the CG method for com- 
plex, symmetric, but non-hermitian matrices. Its algorithm 
is shown by iterative relations among x kl r k and p k vectors 
(k — 1,2,3 • • • ) as, 



xk+i =x k + a k p k , 



n+\ =r k -a k Ap k 



Pk+\ = r k +i 



PkPk, 



(B3) 



(B4) 



(B5) 



where a k = r T k r k /p T k Ap k and p\ k = r T k+l r k+l /r T k r k ( Note that 



a k r k r k jp\Ap k and ft ^ r T k+l r k+l / 'r\r k ). Initial conditions 
are Oo = 1, ft = 0, xq = and ro = b. As iteration number 
k increases, the norm \r k \ of residual vector r k decreases. The 
convergence criterion is given for |rjt|/|fe|. If this convergence 



10 



condition is fulfilled, we can obtain numerically approximated 
solution x. 

Next we consider a series of shifted linear equations as 

(A - al)x a = b, (B6) 

where a is a complex number and 7 is a unit matrix. If we 
start above iteration from xq = 0, the k-th residual vector r° 
of the COCG method for Eq. ( IB6b can be proven to be propor- 
tional to the k-th residual vector r# of the COCG method [7] 
for Eq. (EB (i-e., Eq. dB6| with a = 0); 

r k=^ k , (B7) 

where n? is a proportional coefficient and satisfies following 
iterative relations as, 

< +1 = (l + n^ + ^l^-it,), (B8) 
< = ^S-Uk, (B9) 



These iterative relations can be derived [7] from an invari- 
ance property of two Rrylov subspaces concerning Eqs. dBll ) 
and (IB6b . The former Krylov subspace is generated by the 
iteration of the CG method, that is, 



span{b,Ab,A 2 b,---}. (Bll) 



By shifting A as A — al, the latter Krylov subspace becomes, 



span{b, (A - (?I)b, (A - (Jl) 2 b, ■■■}. (B 12) 



This subspace is the same as that defined in (IB lib . 
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